Crystallography on Curved Surfaces 
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We study static and dynamical properties that distinguish two dimensional crystals constrained 
to lie on a curved substrate from their flat space counterparts. A generic mechanism of dislocation 
unbinding in the presence of varying Gaussian curvature is presented in the context of a model 
surface amenable to full analytical treatment. We find that glide diffusion of isolated dislocations 
is suppressed by a binding potential of purely geometrical origin. Finally, the energetics and biased 
diffusion dynamics of point defects such as vacancies and interstitials is explained in terms of their 
geometric potential. 
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I. INTRODUCTION 

The physics of two dimensional crystals on curved sub- 
strates is emerging as an intriguing route to the engineer- 
ing of self assembled systems such as the "colloidosome" , 
a colloidal armor used for drug delivery [l|, or devices 
based on ordered arrays of block copolymers which are 
a promising tool for "soft lithography" [3, Q. Curved 
crystalline order also affects the mechanical properties of 
biological structures like clathrin-coated pits [J, |3, |g or 
HIV viral capsids 0, 1^ whose irregular shapes appear 
to induce a non-uniform distribution of disclinations in 
their sheU [9|. 

In this article, we present a theoretical and numerical 
study of point-like defects in a "soft" crystalline mono- 
layer grown on a rigid substrate of varying Gaussian 
curvature. We suppose that the monolayer has lattice 
constant of order, say, 10 nm or more. The substrate can 
then be assumed smooth, as would be the case for mono- 
layers composed of di-block copolymers 0, Q. Discli- 
nations and dislocations are important topological de- 
fects that induce long range disruptions of orientational 
or translational order respectively jlO. 11] . Disclinations 
are points of local 5- and 7-fold symmetry in a triangular 
lattice (labeled by topological charges g = ±^ respec- 
tively), while dislocations are disclination dipoles char- 
acterized by a Burger's vector, b, defined as the amount 
by which a circuit drawn around the dislocation fails to 
close (see Figure ^ inset). Other point defects such as 
vacancies, interstitials or impurity atoms create shorter 
range disturbances that introduce only local stretching 
or compression in the lattice (see Figure ISJ. Such de- 
fects are important for particle diffusion and relaxation 
of concentration fluctuations. 

These particle-like objects interact not only with each 
other, but also with the curvature of the substrate via a 
one-body geometric pot ential that depends on the par- 
ticular type of defect |l3, [i3l ■ These geometric potentials 



are in general non — local functions of the Gaussian cur- 
vature that we determine explicitly here for a model sur- 
face shaped as a "Gaussian bump" . An isolated bump of 
this kind models long wavelength undulations of a litho- 
graphic substrate, has regions of both positive and nega- 
tive curvature, and yet is simple enough to allow straight- 
forward analytic and numerical calculations. The pres- 
ence of these geometric potentials triggers defect unbind- 
ing instabilities in the ground state of the curved space 
crystal, even if no topological constraints on the net num- 
ber of defects exist. Geometric potentials also control the 
dynamics of isolated dislocations by suppressing motion 
in the glide direction. Similar mechanisms influence the 
equilibrium distribution and dynamics of vacancies, in- 
terstitials and impurity atoms. 



II. BASIC FORMALISM 

The in-plane elastic free energy of a crystal embed- 
ded in a gently curved frozen substrate, given in the 
Monge form by r(x, y) = (a;, y, h{x, y)), can be expressed 
in terms of the Lame coefficients /i and A 14] 

F^\ jdA ay (x)uy (x) = j dA U u^ (x) + ^ u2,,(x) j , 

(1) 
where x = {x, y} represents a set of standard carte- 
sian coordinates in the plane and dA = dxdy^/g, 

where ^ = y'l + (d^h)^ + (dyh)^. In Eq.© the 



stress tensor o-y (x) 



,(x) 



written in terms of the strain tensor 
^[diUj{'x.) + djUi{x.) + Aij{x.)]. The latter contains an 
additional term Aij(x) = dih{x)djh{x) (compared to its 
flat space counterpart) that couples the gradient of the 
displacement field Mi(x) to the geometry of the substrate 
as embedded in the gradient of the substrate height func- 
tion /i(x). We will often illustrate our results for a single 
Gaussian bump whose height function is given by 



hi^) 



axge 



(2) 
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where r = — is the dimensionless radial coordinate [la . 
The aspect ratio a provides a dimensionless control pa- 
rameter to measure deviations from flatness. 



In Eq.Q and the rest of this paper, we adopt the ordi- 
nary flat space metric (e.g. setting dA « dxdy) and ab- 
sorb all the complications associated with the curved sub- 
strate into the tensor field Aij (x) that resembles the more 
familiar vector potential of clcctromagnetism. Like its 
electromagnetic analog, the curl of the tensor field Aij (x) 
has a clear physical meaning and is equal to the Gaussian 
curvature of the surface G'(x) = —eiiejkdidkAij{'x.) where 
en is the antisymmetric unit tensor {cxy = —^yx — 1) 
[13. The consistency of our perturbative formalism to 
leading order in a is demonstrated in Appendix ^ (Sup- 
plemental Materials). 

Minimization of the free energy in Eq. JQ) with respect 
to the displacements Mi(x) naturally leads to the force- 
balance equation, diaijlpc) = 0. If we write crij(x) in 
terms of the Airy stress function x(x) 



,(x) = euejkdidkxi'^) 



(3) 



then the force balance equation is automatically satisfied, 



diCTij = ejfc9fc[9i,52]x(x) = 



(4) 



since the commutator of partial derivatives is zero. How- 
ever, we must be able to extract from x(^) the cor- 
rect Uy(x) that incorporates the geometric constraint 
of Eq.lO and accounts for the presence of any defects. 
These requirements are enforced by solving a bi-harmonic 
equation for x(x) whose source is controlled by the vary- 
ing Gaussian curvature of the surface G(x) and the defect 
density 5'(x) [Hill 



-A2x(x) = 5(x) 



G(x) 



(5) 



In the special case S'(x) = 0, we denote the solution of 
Eq.© as x'^(x) where the superscript G indicates that 
the Airy function and the corresponding stress tensor 
CTj^ (x) describe elastic deformations caused by the Gaus- 
sian curvature G(x) only, without contribution from the 
defects. The stress of geometric frustration, a^{x), is a 
non — local function of the curvature of the substrate and 
plays a central role in our treatment of curved space crys- 
tallography. The Young modulus Y = ^10^ 113, |l3 
naturally arises upon recasting (up to boundary terms) 
the free energy in Eq.lQJ as a simple functional of the 
scalar field x(x) 



F^^ IdAiAxi^))' 



(6) 



The source S'(x) for a distribution of iV unbound discli- 
nations with "topological charges" {qa = i^} and M 
dislocations with Burger's vectors {b*^} reads [lO| 



5(x) 



N 

E 

Q = l 



fc(5(x,x") + ^e,,&f5,5(x,x'5) 

/3=1 



(7) 



where |b| is equal to the lattice constant a for disloca- 
tions with the smallest Burger's vector. For A'' isotropic 
vacancies, interstitials or impurities, we have |l(1 | 



1 ^ 
5(x) = -^r!„A<5(x,x") 



(8) 



where iln 



is the local area change caused by includ- 



ing the point defect at position x" . 

It is convenient to introduce an auxiliary function V^ (x) 
that satisfies the Poisson equation AV{x.) = G(x) and 
vanishes at infinity where the surface fiattens out. In 
order to determine the geometric potential, C(x"), of a 
defect at x" we integrate by parts twice in Eq. 10) and use 
Eq.© to obtain A^x(x) and x(x) in terms of the Green's 
function of the biharmonic operator. The geometric po- 
tential (on a deformed plane flat at inflnity) follows from 
integrating by parts the cross terms involving the source 
and the Gaussian curvature with the result 



C(x") = -Y fdA'S{x') fdA-^ T/(x) 



(9) 



This formula ignores defect self-energies |31| and needs to 
be supplemented by boundary corrections, as discussed 
in Appendix 151 f Supplemental Materials). 

In the following sections, we explicitly show that our 
results derived from Eq. (|5J) can also be obtained by means 
of more heuristic arguments. According to this point 
of view, defects introduced in the curved space crystal 
are local probes of the preexisting stress of geometric 
frustration cr*^(x) to which they are coupled by intuitive 
physical mechanisms such as Peach-Koehler forces. 



III. 



GEOMETRIC FRUSTRATION 



We start by calculating the energy of a relaxed defect- 
free two dimensional crystal on a quenched topography. 
In analogy with the bending of thin plates we expect 
some stretching to arise as an unavoidable consequence 
of the geometric constraints associated with the Gaus- 
sian curvature 14] . The resulting energy of geometric 
frustration, Fo, can be estimated with the aid of Eq.© 
which, when S'(x) — 0, reduces to a Poisson equation 
whose source is given by T^(x) 



Ax°(x) = -y(x)+if;,,(x) 



(10) 



where i?fl(x) is an harmonic function of x parameterized 
by the radius of the circular boundary R. As discussed in 
Appendix^ iJj^(x) vanishes in the limit R ^ xq ii free 
boundary conditions are chosen. Upon using the general 
definition of the Airy function (Eq. UJ, we obtain 



a^ (x) = Ax^(x) = -ry(x) . 



(11) 



For a surface with azimuthal symmetry, like the bump, 
the only non vanishing components of the stress tensor 



of geometric frustration a^^ (x) read 



.M 



ag{r) 



2^,G 



1 d^x 



1 dx^ 
Xq r dr 



(12) 



(13) 



and Poisson's equation can be readily solved upon apply- 
ing Gauss' theorem with the Gaussian curvature G{r) ~ 



5 as a source [li 



Vir) 



4 1 iri rfrV'G(/') = -i 



2 -r^ 

a e 



(14) 

The extra factors of xq in the previous equations arise 
from expressing our results in terms of the dimensionless 
radial distance r. 

Upon substituting Eg. lfTUI) in Eq.®, we have 



Y 

"2" 



Foc^- IdA V{rf 



2^,4 



YnxQa 
64 



(15) 



The result in Ea. (|15|) is valid to leading order in a, con- 
sistent with the assumptions of our formalism after tak- 
ing the limit R ^ xq (see Appendix B for finite size 
corrections in small systems). For an harmonic lattice, 
Y ~ -j=k, where k is an effective spring constant that can 
be extracted from more realistic inter-particle potentials 
|23j . For colloidal particles, k is typically two orders of 
magnitude larger than kBT/aF', where T is room tem- 
perature [l9j . Our numerical calculations of Fq in fixed- 
connectivity harmonic solids are in good agreement with 
the small a expansion in Ea. (|15|l as long as the aspect 
ratio a is around 1/2 or lower (see Appendix 1X1 for a nu- 
merical plot of Fq versus a) . An immediate implication 
of the geometric frustration embodied in Eq. (|14() and 
(|15|1 is that nucleation of crystal domains on the bump 
will take place preferentially away from the top in regions 
where the surface flattens out. 



IV. GEOMETRIC POTENTIAL FOR 
DISLOCATIONS 

The energy of a two dimensional curved crystal with 
defects will include the frustration energy, the inter- 
defect interactions (to leading order these are unchanged 
from their flat space form, see Appendix , possible 
core energies and a characteristic, one-body potential of 
purely geometrical origin that describes the coupling of 
the defects to the curvature given by Eq.l^. The geomet- 
ric potential of an isolated dislocation, C(x) = D{x,6), 
is a function of its position as well as of the angle 9 that 
the Burger's vector b forms with respect to the radial 
direction (in the tangent plane of the surface. See Figure 
n) Upon setting all q^ — Q in Eq.Q and substituting 



into Eq.®, we obtain, for an isolated dislocation, the 

M 



resulting function D{r = — , ( 
Dir,9) = -Yb,e,,d,JdA'^ (v{^') 



Y bxQ —- sin ( 

8 



-1 



? 2 
4i?2 



|) ^[10) 



In view of the azimuthal symmetry of the surface. Gauss' 
theorem as expressed in Ea. ll4|l . was employed in deriv- 
ing the second equality in Ea. (|16|l which is a function 
only of the dimensionless radial coordinate r. The first 
term in Ea. H16(l corresponds to the R ~* 00 geometric 
potential obtained from Eq. @ , while the second term is 
a finite size correction arising from a circular boundary 
of radius R (see Appendix IBl for a detailed derivation). 
Ea. p6|l is valid to leading order in perturbation theory, 
consistent with the small a approximation adopted in 
this work (see Appendix IX|) . 

In Fig. ^we present a detailed comparison between the 
theoretical predictions for the geometric potential 



D(r,e) 
Ybxn 



plotted versus r = — as continuous lines and numeri- 

^Q . . . . 

cal data from constrained minimization of an harmonic 
solid on a bump with a = 0.5, under conditions such that 
R :^ Xo :$> a. (See Appendix |0 for a discussion of our 
numerical approach). The lower and upper branches of 
the graph are obtained from Eq. ((TC|l by setting 9 = ±^ 
and letting -^ equal to 4 (blue curve) or 8 (red curve). 
Indeed, the (scaled) data from simulations with different 
choices of — (see caption) collapse on the two master- 
curves according to their ratio of ^ . Two different curves 
arise because the dislocation interacts with the curvature 
directly and via its image. The image-mediated interac- 
tion is given by the R dependent term in Ea. ljKil) . The 
sin 6* dependence of D(r, 9) on the direction of the Burg- 
ers vector is revealed by Fig. ^ since the upper branch 
of the graph corresponding to the unstable equilibrium 
9 ^ —^ is approximately symmetric to the lower one 
corresponding to 9 = ^ . 

The analogy between the geometrical potential of the 
dislocation and the more familiar interaction of an elec- 
tric dipole in an external field can be elucidated by re- 
garding the dislocation as a charge neutral pair of discli- 
nations whose dipole moment qdi = tij bj is a lattice vec- 
tor perpendicular to b that connects the two points of 
5 and 7- fold symmetry. The geometric potential, U{r), 
of a disclination of topological charge q interacting with 
the Gaussian curvature satisfies the Poisson equation 
AU{r) = —qV{r) as can be seen by substituting the 
source in Eq.^ with all b'^ = into Eq.(|5J . For small r, 
positive (negative) disclinations are attracted (repelled) 
from the center of the bump by the integrated back- 
ground source V{r) which increases for r < 1 like a^r^ 
and is multiplied by the 2D electric field -, resulting in 
a geometric force that increases linearly in r. If the posi- 
tive disclination within the dipole is closer to the top, the 
force it experiences will be opposite and slightly less than 
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FIG. 1: The dislocation potential D{r,0 = ±f) in Ea. l(T()jl 
(including boundary corrections) is plotted as a continuous 
line for a Gaussian bump parameterized by a = 0.5 in the 
limit R » Xo » a. Open symbols represent the numeri- 
cal minimization of a fixed connectivity harmonic model for 
which the separation of the +/- disclination pair represent- 
ing the dislocation is fixed while varying the dimensionless 
radial distance — of the dislocation from the bump. An an- 
gle 9 = ^ indicates that the + disclination is closer to the 
top of the bump. Lower branch: 9 = tt/2, R/xq =4 (blue), 

8 (red); and xo/a = 10 (Q), 20 (A), 40 (D). Upper branch: 

9 — —n/2, xq/R = 8, xo/a = 10. Inset: A schematic view 
of a dislocation with filled symbols representing six- (circles), 
five- (diamond) and seven- (square) coordinate particles. Also 
depicted are the two rows of extra atoms emanating from the 
five-coordinated particle. The Burger's vector is shown as 
a red arrow, completing the circuit around the dislocation 
(dashed line). 



the one experienced by the negative disclination that is 
further away from the top. As a result a "tidal" force will 
push the dislocation (as a whole) down hill as shown by 
the lower branch of the plot of Fig. ^ For large r, how- 
ever, the source V{r) saturates and the attractive force 
exerted on the positive disclination wins and drags the 
dislocation towards the bump. The minimum of the ge- 
ometric potential occurs at |xmm| ~ 1-1 xq, close to the 
circle of zero Gaussian curvature |x| — xo, as a result of 
the competing interactions of the two disclinations com- 
prising the dislocation. If the orientation of the discli- 
nation dipole is flipped, the geometric potential changes 
sign. Similarly, if the sign of the Gaussian curvature is 
reversed, the sign of the geometric interaction flips. As a 
result, a dislocation close to a saddle will have its Burger 
vector oriented so that the closest disclination to the sad- 
dle is 7-fold coordinated. 

The physical origin of the dislocation potential can 
be understood heuristically by considering dislocation 
climb, or the motion of the dislocation in the radial direc- 
tion (Figure [2J|. According to standard elasticity theory, 
a dislocation in an external stress field aij (x) experiences 




FIG. 2: Schematic of dislocation climb for the case of a square 
lattice, (a) Vacancy-mediated dislocation climb depicted in 
four steps: (i) Initial configuration consisting of an isolated 
dislocation (grey) (bonds in the third row are highlighted with 
dashed lines) ; (ii-iii) Motion of a vacancy (white) to the dis- 
location; (iv) Motion of the dislocation in the upward radial 
direction to create a configuration analogous to (i). Compari- 
son of the third row of atoms in panel (i) and (iv) reveals that 
this process results in over-stretching of the bonds by a length 
of order of the lattice spacing a. (b) On a bump (outline of 
patch shown), the process represented in (a) is energetically 
unfavorable close to the top of the bump because there is a 
tensile stress a^^ (force per unit length indicated by arrows) 
that hinders further stretching of the bonds and tries to push 
the dislocation downward with a Peach-Koehler force a^^b 
where the Burger's vector |b| ft; a. The force changes sign 
at a distance of approximately xq from the top of the bump 
(black ring) where a^^ switches sign or if the direction of b 
is reversed. 



a Peach-Koehler force, f(x), given by //c(x) = efcj6icry(x) 
|16l llTj . Similarly, a dislocation introduced into the 
curved 2D crystal will experience a Peach-Koehler force 
as a result of the pre-existing stress field of geometric 
frustration a^, (x" ) whose non-diagonal components van- 
ish. This interpretation is consistent with the geometric 
potential derived in Eq. H16|) , provided we use Eq. H1U|I to 
write -D(x) — hieijdjx'^ ■ With b along its minimum ori- 
entation (azimuthal counter-clockwise), we obtain a ra- 
dial Peach-Koehler force of magnitude /(r) = —ba9Jr) 

.aD(r)-_A^!x^(seeEq.C3). 



that matches 



Xq dr 



dr'^ 



DISLOCATION UNBINDING 



VI. GLIDE SUPPRESSION 



If the 2D crystal is grown on a substrate which is suffi- 
ciently deformed, the resulting elastic strain can be par- 
tially relaxed by introducing unbound dislocations into 
the ground state [Tj, [13 • Here we present a simple esti- 
mate of the threshold aspect ratio, ac, necessary to trig- 
ger this instability. Boundary effects will be ignored in 
what follows by letting i? ^ cx) in Ea. H16|l . 

Consider two dislocations located at xi and X2 a dis- 
tance of approximately xq from the center of the bump 
on opposite sides (see the inset of Fig. O. Their 
disclination-dipole moments are opposite to each other 
and aligned in the radial direction so that the two (an- 
tiparallel) Burger's vectors are perpendicular to the sep- 
aration vector in the plane X12 = Xi — X2. In this 
case, the interaction between the dislocations reduces to 

V12 ~ £ ^^ In(^) Eini- The instability occurs 
when the energy gain from placing each dislocation in 
the minima of the potential D(x.) given by Ea. (|16|l out- 
weights the sum of the work needed to tear them apart 
plus the core energies 2Ec- The critical aspect ratio at 
the threshold which results is given by 



i-(i^) 



(17) 



where b' = | e~ ^^'^ and c « 1/2. Note that the core 
energy E^ is determined by the microscopic physics of 
the particular system under study. 

In Fig. |2t we present a comparison between Eq. p7|l 
and numerical results. For each value of ^, the corre- 
sponding ttc is obtained numerically by comparing the 
energy of a lattice without defects to the configuration 
with the two dislocations in their equilibrium positions. 
This interpretation for the origin of the instability is cor- 
roborated by the (numerical) strain energy density plots 
of Fig. O) andO; , where it is shown that introducing the 
pair of dislocations reduces the strain energy density on 
top of the bump at the price of creating some large, but 
localized strains around the dislocation cores where m^ 
diverges. In the continuum limit & <C xq, very small de- 
formations are enough to trigger the instability. This is 
the regime in which our perturbation treatment applies. 
As a is increased even further, a cascade of dislocation- 
unbinding transitions occurs involving larger numbers of 
dislocations and more complicated equilibrium arrange- 
ments of zero net Burger vector. For sufficiently large 
aspect ratios, we expect that the dislocations tend to 
line up in grain boundary scars similar to the ones ob- 
served in spherical crystals [l| . This scenario is consistent 
with preliminary results from Monte Carlo simulations in 
which the fixed-connectivity constraint is lifted and more 
complicated surface morphologies are considered [l8|. 



The dynamics of dislocations proceeds by means of two 
distinct processes: glide and climb. Glide describes mo- 
tion along the direction defined by the Burger's vector; in 
flat space glide requires a very low activation energy and 
is the dominant form of motion at low temperature (see 
Fig. 0^ inset). Climb, or motion perpendicular to the 
Burger vector requires diffusion of vacancies and intersti- 
tials (Figure [SJ and is usually frozen out relative to glide 
which involves only local rearrangements of atoms. On a 
curved surface, the geometric potential D{r, 9) imposes 
constraints on the glide dynamics of isolated dislocations, 
in sharp contrast to flat space where only small energy 
barriers are present due to the periodic Peierls potential 

M 

As the dislocation represented in the inset of Fig. 0^ 
moves in the glide direction, it experiences a restoring po- 
tential generated by the variation of the (scaled) radial 
distance and the deviation from the radial alignment of 
the dislocation dipole. For a small transverse displace- 
ment, y, the harmonic potential U{y) = 2^dV^ is con- 
trolled by a radial, position-dependent effective spring 
constant, kd which depends on the radial coordinate r. 
Upon expanding Ea. H16|l to leading order in y and a, we 
obtain 



kdir) 



a'^Yb (1- (l + r2)e- 
4 xo I r^ 



(18) 



The harmonic potential 2^d.y^ is shown in Fig. 2K where 
the data obtained from numerical minimization of the 
harmonic lattice is explicitly compared to the prediction 
of Ea. H18|) for different values of the aspect ratio. Note 
that the effective spring constant plotted in Fig. ^ van- 
ishes in the limit — > but can still be important for 

small systems since Yb^ can be of the order of hundreds 
oi ksT = f3~^ (see section ^QJ. The confining potential 
plotted in Fig. ^^ is similar to the one experienced by a 
dislocation bound to a disclination [l9l |2(]| . 

The resulting thermal motion in the glide direction of 
dislocations in this binding harmonic potential can be 
modeled by an over-damped Langevin equation for the 

where 



-2iik^t 



Pka 



glide coordinate y, leading to (Ay ) = 
/i is the dislocation mobility [l^ ^^. In the case of a 
bump geometry, the effective spring constant kd can be 
evaluated using Ea. (fTS|l . We emphasize, however, that 
the glide suppression mechanism considered here is not 
caused by the interaction of the dislocation with other 
defects but purely by the geometric interaction with the 
curvature of the substrate. 



VII. VACANCIES, INTERSTITIALS AND 
IMPURITIES 

We now turn to a derivation of the geometric potential, 
/(x), for interstitials, isotropic vacancies and impurities. 




FIG. 3: (a) Dislocation unbinding critical aspect ratio, a^, as a function of ^. The theoretical estimate (1171 is plotted vs. xq 
for core energies Ec — (dashed) and the best fit Ec = O.lYb^ (solid). Numerical values (circles) were obtained as described in 
SectionEland Appendix^ Inset: particle configuration projected in the plane for a dislocation pair straddling the bump. The 
extra atoms associated with the dislocations are highlighted with black lines, (b) Log of the numerical strain energy density on 
a bump for the configuration shown in the inset, (c) The same quantity for a defect-free bump. Both plots were constructed 
numerically with xq = 10, a = 0.7 > Uc- Red represents high strain. 




(a) 




FIG. 4; Dislocation glide, (a) Filled circles represent numerical glide energies vs. distance along the glide direction y (units of 
a) for xo/a = 10, R/xq = 8, r = 0.5 (at y = 0): a — 0.1 (dark blue), 0.3 (blue), 0.5 (orange) and 0.7 (red). Solid lines represent 
E{y) = ^kdy^, with kd determined from Eq. HSU . The energy is scaled by 10~*. Inset: schematic of a dislocation (red) on a 
Gaussian bump. The glide path is highlighted in orange, (b) The curvature- induced glide-suppression spring constant kd{r) 
1181 1 (solid line) is plotted versus scaled in-plane distance from the top of the bump r — |x|/a:o for xo = 10a, a = 0.5. Open 
symbols represent numerical results found by fitting similar data as in (a) to parabolic curves in the glide coordinate. The 
ordinate axis is scaled by lO"'^. 



Inspection of Fig. [Slreveals that an interstitial (vacancy) 
can be viewed either as the product of locally adding (re- 
moving) an atom to the lattice or as a composite object 
made up of three disclination dipoles. In order to derive 
/(x), we substitute the source term of Eq.® into Eq.Q 



and integrate by parts twice. The result reads 



(19) 



where boundary terms and a position independent nucle- 
ation energy have been dropped. The constant 57 repre- 
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FIG. 5: The scaled potential energy I{r)/Ya? of interstitials 
(bottom) and vacancies (top) is plotted versus the scaled dis- 
tance r for the same bump as Fig. Solid lines are obtained 
from Eq. ()19^ . where the area changes of interstitials and 
vacancies Q.f and fl" were fit to the data. Sample configu- 
rations are plotted in the insets, where diamonds, circles and 
squares represent five, six and seven fold coordinated parti- 
cles respectively. Filled and unfilled circles in the interstitial 
plot represent two distinct but energetically equivalent orien- 
tations of the disclination dipoles comprising the interstitial. 
They differ by swapping the 5 and 7- fold disclinations, which 
rotates the orientation of the filled, triangular plaquette by 



sents the area excess or deficit associated with the defect. 
In Fig. |31 a comparison of Eq. (|19|l with the results ob- 
tained from mapping out /(r) numerically is presented. 
The area changes fi^ and r2„ for interstitials and vacan- 
cies respectively were fit to the numerical data. We find 
that r^i, is negative and greater in magnitude than fi^. 
The large r behavior of I{r) indicates that the core en- 
ergy of a vacancy in flat space is greater than the one 
of an interstitial for the harmonic lattice. Interstitials 
tend to seek the top of the bump whereas vacancies are 
pushed into the flat space regions. From the deflnition 
of V{r) in Ea. (|14|l . we deduce that an interstitial (va- 
cancy) is attracted (repelled) by regions of positive (neg- 
ative) Gaussian curvature similar to the behavior of an 
electrostatic charge interacting with a background charge 
distribution given by —G{r). The function V(r) controls 
the curvature-defect interaction for other types of defects 
such as disclinations in liquid crystals and vortices in ^He 
films [iajllj. The expression for V{r) in Ea. H14(l reveals 
that J(r) is indeed a non local function of the Gaussian 
curvature determined as in electrostatics by the appli- 
cation of Gauss' law. Thus, the vacancy potential on 
a bump does not reach a minimum at the point where 
the Gaussian curvature is maximally negative but rather 
at infinity where the integrated Gaussian curvature van- 
ishes. 



We now argue heuristically how a localized point de- 
fect can couple non-locally to the curvature and in the 
process we provide an alternative justification of Eq. 
(|19|1 analogous to the informal derivation of the dislo- 
cation potential. The energy cost of a defect at x", 
/(x"), due to local compression or stretching in the 
presence of an arbitrary elastic stress tensor cry(x) is 
given by /(x") = p{-x.'^)5V where 5V is the local vol- 
ume change and p(x") is the local pressure related to 
the stress tensor via ^^(x") = —p{x")Sij for the case 
of an isotropic stress[Ij|. In two dimensions, we have 
/(x") = -crfcfc(x")^. We recover the result in Ea. (P|l . 
by assuming that the local deformation Q" (induced by 
the nucleation of a point defect) couples to the preexist- 
ing stress of geometric frustration a^^. = —YV{x.) (see 
Ea. (|ll|l ). which is a non-local function of the Gaussian 
curvature. Elastic deformations created by the geometric 
constraint throughout the curved two dimensional solid 
are propagated to the position of the point defect by force 
chains spanning the entire system. The point defect can 
then be viewed as a local probe of the stress field that 
does not measure the additional stresses induced by its 
own presence. 

Note that the geometrical potential of an isotropic 
point defect is unchanged if we swap the 5 and 7-fold 
disclinations comprising it (corresponding to a rotation of 
the point defect by ^ around its center) , as demonstrated 
for interstitials in the lower branch of Fig. [S] Contrast 
this situation with the clear dipolar character of the dis- 
location potential plotted in Fig. Q] The more compli- 
cated case of non-isotropic point defects appears to still 
be captured qualitatively by Ea. ljiyil : this was explicitly 
checked by plotting the geometric potential of "crushed 
vacancies" |2^ . which have both a lower symmetry and 
a lower energy than their isotropic counterparts. 

An arbitrary configuration of weakly interacting point 
defects will relax to its equilibrium distribution by diffu- 
sive motion in a force field f (r)= — V/(r). In overdamped 
situations, this geometric force leads to a biased diffusion 
dynamics with drift velocity v ~ (3\f\a— ^ — —, where 
/3 = -j-ijT and D is the defect diffusivity [la|- Eventually, 
a dilute gas of point defects equilibrates to a non-unform 
spatial density proportional to e~^-^^'^\ 
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The supplemental materials are organized in the form 
of three Appendices A, B, and C. 



APPENDIX A: PERTURBATION THEORY OF 
CURVED CRYSTALS 

The starting point of our perturbative analysis of 
curved crystalline order is Eq.Q which incorporates 
Gaussian curvature by adding an extra source to the de- 
fect term but it is nonetheless written using a flat space 
metric. To underscore the subtleties involved we write 
a covariant generalization of the force balance equation, 

C).fTy(x)=0, in 



;^|T(V?^])-r>l = o. 



(Al) 



where g is the determinant of the metric tensor gij and 
r^j is the Christoffel symbol. Ea. ljAljl can be concisely 
written as ct'-.j = 0, where the semicolons indicate the 
covariant derivatives Di |2J|. 

Strictly speaking, this equation cannot be solved by 
using the flat space trick of writing the stress tensor in 
terms of the Airy function because of a distinctive prop- 
erty of curved space: torsion |25Ll2q. An arbitrary vector 
T parallel transported around a closed loop on a surface 
of non vanishing Gaussian curvature is rotated from its 
original orientation |32l |. In differential form, this con- 
straint reads 



[A,i?,]T, = G(x)7,,7rr„ 



(A2) 



where 7^^ = -j= denotes the antisymmetric tensor on a 

curved surface [26]. Upon substituting Tm = 7m" dn x(x) 
into Ea. (|A2|) . we see that the covariant analogue of 
Eq. Q) does not hold because the commutator of covari- 
ant derivatives (known as the torsion tensor) does not 
vanish. 

We can nonetheless make progress by noting first that 
the Gaussian curvature of a bump in Eq. I)A2(I is propor- 
tional to % and reads [la 



G{r) = 



c?e-- 



1 



fl -|- cP'T^e ^^ 



(A3) 



As a result, Eq. Q is in fact approximately correct 
for small a. We first consider the case ^(x) — 0, for 
which Eq. Q reduces to one of the two celebrated Foppl- 
von Karman equations describing slightly deformed thin 
plates [ij|. For a frozen substrate, the second Foppl- 
von Karman equation (arising from the variation of the 
normal coordinate h(x,y)) determines the adhesion pres- 
sure (normal to the surface) which is needed to constrain 
the 2D solid to the curved substrate [lj|. The Foppl- 
von Karman analysis rests on a consistent perturbation 
theory in a and predicts that x(x) is O(a^), as can be 
seen from Eq.jS]) with S'(x) — 0. To this order, the com- 
mutator in Ea. (|A2|) indeed vanishes (the leading order 
corrections are 0{a^)) and one is justified in expressing 
crij(x) in terms of an Airy function x(x)- 
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FIG. 6: The geometric frustration energy Fo{a) (Eq. 1151 ) is 
plotted versus a for xo/R = 10/80 (Q) and 10/40 (D). Solid 
{xo/R = 10/80) and dashed (xp /R = 10/40) lines indicate 
plots of Eq. I15II using y = -2=fc |23l . including the finite size 
corrections discussed in Appedix |B] for these two conditions, 
respectively. 



The situation is more delicate in the presence of defects 
because x(^) ^^o longer vanishes as a ^ but tends in- 
stead to the flat space form (see Ref. [l^, UM fo^ explicit 
results). Dimensional analysis reveals that the commu- 
tator in Ea. (|A2|l is now proportional to ^^— for dis- 

interstitials and 



locations and to 



xo xl 



for vacancies, 



impurities where corrections of order In (■^j are ignored. 
Hence, the commutator in Ea. (|A2|l . set to zero in the 
derivation of Eq. Q , appears to be of the same order in 
a as the curvature corrections to x(^) that we wish to 
calculate. The commutator is still small, however, in the 
continuum limit a:o ^ a. The use of Eq.© to study iso- 
lated disdinations to leading order in a is more difficult 
to justify because in this case the commutator can be as 
large as ^-^-^ . However, such an investigation on a sur- 
face with the topology of the plane is of limited interest, 
in view of the large energy cost of isolated disdinations 
(quadratic in R). In this work, we concentrate primar- 
ily on the physics of dislocations, vacancies, interstitials 
and impurities, which (as confirmed by our simulations) 
is adequately described by the formalism embodied in 
Equations 10) and ®. Our analytical results are com- 
pared whenever possible to numerical minimizations of 
an harmonic lattice model; these computations corrobo- 
rate our qualitative conclusions even beyond the limits 
a ^ 1 and a ^ xq for which Equations 10) and (@ are 
strictly valid. For example, Fig. [Sjillustrates how our the- 
oretical predictions of the frustration energy Fq, based on 
Ea. (|15|) . fit numerical data for increasingly larger values 
of a. Similarly, rather good agreement between theoret- 
ical predictions and numerics is obtained for the disloca- 
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FIG. 7: The dislocation potential D(r, 7r/2) in Eq.lQlJ is plot- 
ted as a continuous line for a Gaussian bump parameterized 
by a = 1. Open symbols represent the numerical minimiza- 
tion of a fixed connectivity harmonic model for which the 
separation of the two 5-7 disclinations is fixed while sliding 
the dislocation as a whole radially with respect to the bump: 
e = 7r/2, R/xo =4 (blue), 8 (red); xo/a = 10 (Q), 15 (□), 
20 (A), 30 (0). The R-dependent finite size corrections of 
Appendix f51 are included. 



tion potential D{r, 9) even if a 

m 



1, as illustrated by Fig. 



explicitly determine Hji by integrating 



1 



-Ax''{r) = ~V{r)+Hn 



over the area of the circular disk, with the result 
R rR/^o 



,9x 
dr 



= Yxl / [Hr - V{r)] rdr 

"'0 



(B2) 



(B3) 



Upon substituting Equations (jBlj) and (fn|) into Ea. l|B3|) . 
we obtain by integration 



H, 



4R2" 



(B4) 



which insures that the forces on the boundary vanish. 



Note that cr£c 



1 1 d'x'' 



as a consequence of az- 



imuthal symmetry, so the second boundary condition is 
automatically satisfied. For large system sizes i? ^ xq, 

2 2 

Upon using the general definition of the Airy function 
in the same limit, we obtain 



a,^,(r)=.Ax«(r)^-ry(r)-"'^^° 



4i?2 



(B5) 



Substitution of Eg. ljBSp in Eq.© gives an estimate of 
the stretching energy, Fq, of the defect-free crystal that 
accounts for finite size effects, namely 



APPENDIX B: 



CURVATURE INDUCED FINITE 
SIZE EFFECTS 



In this appendix, we discuss how the geometric poten- 
tials derived above for an infinite system are modified by 
the presence of a circular boundary of radius R on which 
free boundary conditions apply 



,iR) = 



1 

xqR 


\dx] 
dr 



= . 



(Bl) 



where n^ is the unit normal to the circumference of the 
system. We first determine the Airy function x'~^{^) tti^t 
describes elastic deformations caused by the Gaussian 
curvature G{r) only, without any contributions from the 
defects. As explained in Section lllll once x'^i^) is known 
the energy of geometric frustration can be easily calcu- 
lated. We start by fixing the harmonic function Hii{x.) 
introduced in Ea. l|10() upon using the boundary condi- 
tion of Eg.ljBljl. Note that by azimuthal symmetry, the 
Gaussian curvature G{r) and hence Hii{r) are constant 
on the circular boundary. This allows us to conclude that 
the harmonic function is constant everywhere; we denote 
this constant by Hr. The subscript indicates that the 
constant Hji is a function of the system size R. We can 



Fn 



Y 

"2" 



dA 



V{r) 



a Xq 



4i?2 



(B6) 



The graph in Fig. is a numerical plot of Fq versus a 
for — equal to 8 and 4 that corroborates the results of 

xn ^ 

Ea. (jB6|) . This result is obtained from the one quoted 
in the main text by simply performing the substitution 
V{r) -^ V{r) - Hr in Eq.ljTSI). Note that the form of Hr 
was determined by solving for x'^(r) in Ea. HB2|l . Strictly 
speaking, in the presence of defects one should solve the 
full Eq.(|SJl that accounts for the presence of defects by 
means of an extra source term. The solution x(^) would 
not be azimuthally symmetric especially when defects are 
located well away from the center of the bump. However, 
as detailed in Sections IIVI and IVIII the geometric forces 
experienced by the defects at different locations on the 
bump can be easily calculated once cr^{r) is known with- 
out solving the full biharmonic equation provided that 
Xq <C R- The leading finite size correction to a-fjr) was 



indeed calculated in Ea. (|B5l) . 

We can thus make progress in the calculation of the 
defect potentials in the presence of the boundary by sim- 
ply letting V{r) -> V{r) ~ Hr in Eq.®. The defect 
potential now reads 

a^") = -Y JdA'Sii^') JdA^ [Vi^)-HR] (B7) 
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FIG. 8: The dislocation potential D{r, n/2) is plotted versus 
the scaled coordinate r = x/xo, xo = 10a, R = 40a for a = 
(O), 0.05 (A), 0.1 (V) and 0.2 (□). The numerical energy is 
scaled by 10~^a;o, and shifted so that -D(O) = 0. Notice that 
the minimum is gradually lost for a < as ~ 0.1. 



The result for the geometric potential of dislocations 
D{r,9) with finite size effects (see Section Hvj) becomes 



D{r,e) - ^Yhe,,d, JdA' ^ (v{r') 



Y bxo —- sin( 



2 2 

a Xq 
4i?2 



i)%<B«, 



of this work. Here we will only provide some estimates 
of how B(r, 6) (present when R is finite) wins out over 
the geometric interaction D{r,9) in the flat space limit 
a ^ 0. The geometric potential that tries to confine the 
dislocation close to the bmnp, scales like D{r, 6) ^ Ybxo 
(see Ea. (|B8|l V whereas the attraction to the "image dislo- 
cation" is approximately given by B{r,6) « Yb^ [r] ■ 
This choice insures that the boundary interaction is un- 
changed if the position of the defects is flipped x -^ — x 
and suggests that B{r, 9) evaluated at the minimum of 



D{r,6) is given approximately by 



Yb" 



i?2 



(to leading or- 



der in a perturbation expansion in the small parameter 
^). A simple force (or equivalently an energy balance) 
argument provides an estimate of the spinoidal aspect 
ratio Us ~ ^^° below which the minimum of the geo- 
metric potential D{r, 9) is lost to the boundary attraction 
B{r,9). Our estimate agrees with the numerical results 
illustrated in Fig. |H| 



APPENDIX C: NUMERICAL METHODS 

We have complemented our analytical studies with nu- 
merical minimizations of the in-plane stretching energy of 
a triangular lattice of points connected by springs draped 
over the Gaussian bump described in the text. The dis- 
crete stretching energy for a set of points at positions r^ 
is defined by 



F 



discrete 



~ A Z_^ f^fi.iyVfi,!^ 



(Cl) 



fljU 



The second term in D{r, 9) can be viewed as the geomet- 
ric potential of an image dislocation in an infinite system 
(let i? — *■ oo in Ea. H16|l ^ evaluated at the image position 
of |x|' — -^ and with Burger's vector b' = — b. Thus, 
the dislocation interacts with the curvature directly and 
via its image. This choice of an image defect insures that 
the stress normal to the boundary, a^r = -^--g^, van- 
ishes. The curvature induced boundary term included 
in Ea. (|B8|) has no analogue in flat space but its effect is 
evident upon examining the results of our curved space 
simulations presented in Fig. ^for two different choices 
of — equal to 8 and 4 and aspect ratio a = 0.5. The com- 
parison between numerics and the analytical expression 
derived in Eq. IJB8|I confirms the validity of our approxi- 
mation scheme. 

We should stress that the image defect does not guar- 
antee that the second boundary condition Ur^, = is sat- 
isfied when the presence of an off-axis dislocation breaks 
the azimuthal symmetry. However, we argue qualita- 
tively that, as the dislocation approaches the boundary, 
it will be attracted to it as an electrostatic charge would 
to its image. A systematic treatment of this boundary in- 
duced self energy B{r, 9) for arbitrary orientations of the 
Burger vector of the dislocation is outside of the scope 



where r^,^ = |r^ - y^\, r^ = (x^, y^, /i(a;;^,y^)), and a 
is the lattice spacing. The height function h{x,y) de- 
fines the fixed topography of the system and is given by 
h{x,y) — axQC^'^^ "''^ )/2^o. The spring constant matrix 
is fc^,iy — krif^^i,, where the connectivity matrix n^^i, spec- 
ifies that the underlying lattice is triangular, with a fully 
coordinated particle having 6 nearest neighbors (n.n.) 



M;'^ 



1 if fj,,!/ n.n. 
otherwise 



(C2) 



Defects are introduced into the lattice by changing the 
number of nearest neighbors from 6 to 5/7 for +/- discli- 
nations, respectively. Dislocations, interstitials and va- 
cancies are composed of specific configurations of +/- 
disclinations. By coarse-graining Ea. (|Cl|l |23l |. we can 
make an explicit connection to the macroscopic energy 
formula |^ with Young modulus, Y — -y=k and Poisson 
ratio, a — 1/3, from which the Lame coefficients A and 
/i can be determined |l4| . In this work, we use units of 
length and energy such that a = 1 and fc = 1. 

Numerical minimizations of Ea. ljCl|l are performed for 
circular patches of triangular lattice draped over a bump, 
with frozen-in defect configurations. Thus every point in 
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TABLE I: Particles constrained in energy minimizations. 



Configuration 



Number of Particles Fixed 



Particles Fixed 



Particle at (x,y) = (0,0) 
+/- Disclination pair 
+/- Disclination pairs 

Particle at (x,y) = (0,0) 



Defect Free 

Isolated Dislocation 

Two Opposing Dislocations 

Interstitials/Vacancies (I/V) 



Figure^ for example, represents a separate energy mini- 
mization. To perform the minimization, we use the stan- 
dard Fletcher-Reeves (FR) conjugate-gradient method 
|28l |29| in which particles are moved along successive 
directions determined by the gradient of the energy in 



J 



Ea. (|Cl|l with respect to particle coordinates. [33| Af- 
ter taking into account the fact that the z-coordinate is 
determined by h{x,y), the gradient with respect to the 
x-coordinate of particle tj reads 



OF 

dXn 



Mi^ 



M,i> _jiv ^^^^^^ _ ^^^^^^ I ^^^ _ ^^^ ^ [h{x^,,y^j) - h{x„,y^)] - — h{xjj,yrj) 



(C3) 



r 



with 5^,1, the Kronecker delta. To obtain the gradient 
with respect to the y-coordinate, interchange x <-*■ y. 

Convergence in the gradient is achieved when the mag- 
nitude of the gradient of the energy drops below some de- 
fined tolerance, which was set to 10^^ (fc — 1). In order 
to ensure a smooth approach to the energy minimization, 
gradient convergence was accepted only if the | A_F'^'^'^''°*°| 
between the last two iterations of the algorithm was less 



than 10-8 (fc = 1). 

Since there is a non-zero frustration energy for a defect- 
free lattice on a bump, the particles will collectively slide 
off of the bump into flatland if allowed. To prevent this, 
some particles were always fixed during minimization as 
indicated in Table|l| These constraints were implemented 
so that a particle was always located at the center of the 
bump. 
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For dislocations, vacancies and interstitials the position- 
dependent self energies can be ignored compared to C(x°') 
because they are proportional to higher powers of the 
lattice spacing a. 

In fact, the net effect of parallel transport is equivalent to 
monitoring the change of T first in one direction then in 
the other followed by subtracting changes in the reverse 
order. 

The FR algorithm requires the use of a one-dimensional 
minimization algorithm, for which we used the Brent al- 
gorith m [ 2^1 as implemented by the Gnu Scientific Li- 
brary |3(l |. The Brent algorithm requires bounds to be 
placed on the dimensionless parameter A, which we typ- 
ically take to be —5 < A < 10. 



